Analysis of parameter mismatches in the master stability function 
for network synchronization 
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Abstract. - In this letter, we perform a sensitivity analysis on the master stability function 
approach for the synchronization of networks of coupled dynamical systems. More specifically, we 
analyze the linear stability of a nearly synchronized solution for a network of coupled dynamical 
systems, for which the individual dynamics and output functions of each unit are approximately 
identical and the sums of the entries in the rows of the coupling matrix slightly deviate from zero. 
The motivation for this parametric study comes from experimental instances of synchronization 
in human-made or natural settings, where ideal conditions are difficult to observe. 



Introduction. Synchronization of networks of coupled 
dynamical systems has been the subject of intensive re- 
search, see for example the reviews |lHl]- Chaos synchro- 
nization of networked dynamical systems finds applica- 
tions in secure communication [IHT] , system identification 
[SUllj. data assimilation [T21IT3]. sensors [14], information 
encoding and transmission |151I16| , and multiplexing |17j . 
In this framework, the master stability function analysis 
provides a necessary and sufficient condition for the linear 
stability of the synchronous solution. However, most of 
the research on this approach focuses on ideal conditions, 
which are difficult to implement in experiments. 

We consider a typical experimental scenario for a set 
of dynamical systems that are coupled through a network 
to achieve synchronization. We assume that each of the 
elements which constitute the experiment is selected to 
reflect certain nominal characteristics; yet, we allow these 
components to be affected by small mismatches from their 
nominal values. We consider a wide range of possible 
deviations from nominal operating conditions that may 
affect simultaneously the individual units' dynamics, the 
individual units' output functions, and the coupling gains 
among the systems. Another motivation for the proposed 
analysis is the study of the collective behavior of biological 
groups, where individuals are generally different in nature 
and their couplings are typically affected by fluctuations 
about an average or nominal value; see for example '18J. 

We consider the following equations of motion for a set 



of coupled chaotic systems in their nominal conditions 



N 



NOM 



H{xj{t)), i = l,2,...,N, 
(1) 

where Xi G M" is the n-dimensional vector describing the 
state of node i, F : M" — > K" governs the uncoupled dy- 
namics of node i, H : M" — M" is a vectorial output 
function, cr is a scalar gain describing the overall coupling 
strength, and N is the number of nodes in the network. 
The network is deflned by the matrix A^'^^ = {Af^^^^'j, 
describing the coupling from node j to node i. We re- 
fer to equation set ([1]) as nominal, as we assume that it 
corresponds to a given experimental design. A sufficient 
condition for the existence of a synchronized solution. 



Xi(t) = X2{t) 



is that 



E 
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XNit) = Xsit), 



= 0, l = l,...,N, 



(2) 



(3) 



that is, all the row-sumf0 of the matrix A^'-"^'^ are equal 
to zero. In case condition ([3]) is satisfied, a synchronized 



^In what follows, we sometimes refer to the row-sums of a matrix, 
indicating with this terminology the sums of the entries along the 
rows of the matrix. We further comment that the analysis stays 
unaltered if the right hand side of JSJ equals a constant. 
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solution Xs{t) exists that satisfies 

Xs{t) = F{Xs{t)). 



(4) 



We use £^o^'',...,e%'="^' to identify tlie eigenvalues of 
the matrix A^om ^ which are in general complex num- 
bers. Note that © imphes that yl^OA/ j^^^g 

one eigen- 
value, ^^OA/ _ with associated right eigenvector Ijy — 
[1,1,...,!]. 

The linear stability of ([2]) can be assessed by using the 
master stability function |19 1 l20 j . Within this framework, 
the synchronous solution is stable if the maximum Lya- 
punov exponent associated with the parametric equation 



j{t) [DF{xs{t)) + cDH{xS))]l{t) 



(5) 



is negative for every c = crv^^ , k — 2,...,iV, where 7 
is an n-dimensional vector. Then it is possible to asso- 
ciate a master stability function to Eq. ([5]), which yields 
the maximum Lyapunov exponent of ([5]) as a function of 
the parameter c. Thus stability of the synchronized solu- 
tion can be assessed for any given network described by 
Eqs. ([ij by verifying that the corresponding eigenvalues 
a£^'~'^^ , k = 2, N, are within the region of the complex 
plane for which the master stability function is negative. 

Problem statement. The assumptions underlying the 
set of equations ([T]) are that: 

(i) the individual units are all described by the same 
dynamics Xi{t) = F{xi{t)); 

(ii) the systems' outputs are all described by the same 
function H; 

(iii) the sums of the rows of the matrix ^^oa/ ^^^.^ 
zero, that is, condition ([3|) is verified at each node i = 
l,...,iV. 

While assumptions (i)-(iii) can be easily reproduced in 
a numerical simulation, their practical implementation in 
experiments is challenging. Qualitatively good satisfaction 
of (i)-(iii) in experimental instances of synchronization of- 
ten requires fine tuning (2TH28] . In [27H29] . an adaptive 
strategy to dynamically preserve synchronization in the 
presence of slow a-priori-unknown time-variations of the 
couplings is proposed. Though such strategy is able to 
preserve condition (iii) in the presence of external pertur- 
bations, the row-sums of the coupling matrix are typically 
non zero over the time scale of the adaptation. 

In [30,31], assumption (i) is removed and the effect of 
small mismatch of the individual units is considered. That 
is, these works consider the case where F in Eq. (1) is re- 
placed by Fi and the difference between F and Fi is small. 
In this letter, we extend the considerations of [30,31] to 
simultaneously allow for deviations from the exact satis- 
faction of all three of the assumptions (i), (ii), and (iii). 
Namely, we assume that (i),(ii), and (iii) are nominal de- 
sign conditions, which might not be exactly reproduced in 
an experiment. We show that if all the mismatches are 
small as compared to the nominal conditions, the linear 
stability of the nearly synchronized solution can be studied 
by using an extended master stability function. Moreover, 



when the nearly synchronous evolution is stable, the mis- 
matches introduce forcing terms in the parametric equa- 
tion that maintain the network in a state of approximate 
synchronization. 

To take into account approximate, rather than exact 
satisfaction of (i), (ii), and (iii), we rewrite the network 
equations in the form 



N 



Xi(t) = F{xi(t),mi 



+ c7^Ayii-(a;,(t),Pj), (6) 
i=i 



i — 1,2,...,A^, where represents the coupling from 
node j to node i, is a parameter used to identify vari- 
ations of the dynamics at each node i, and pi is a param- 
eter of the output function of each node i. We assume 
that TO; — fh + Srui, where m = rrii and Snii is a 

small mismatch. Similarly, we write Pi = p + 5pi , where 
p = N~^J2iPi ^-nd Spi is a small mismatch. Note that 
by construction Srui — and ^Pi = 0- The ele- 
ments Aij^s represent imperfect realizations of the nom- 
inal couphngs 4^'='*^'s, that is, Aij = Af°^ + 6A,j, 
i,j = 1, N, where SAij is a small mismatch. In general, 
in the presence of deviations of the Aij^s from their nom- 
inal values, it is not possible to write a condition equiva- 
lent to and thus to extend directly the master stability 
function formalism. For small JA^ 's, we can write 



where 



Aij = SAij = Sa + Stti 



Sa = N^^ A,j ^ N^^ 6Ai 



(7) 
(8) 



is the average sum of the rows of the matrix A and. 



Sa,, = { > SAu ] - Sa = 



(^M,,)-iV-i5]M,, (9) 



3 3 V 

is a small deviation. The deviations Sai are calculated 
with respect to the average row-sum Sa, hence they have 
zero sum, that is, Saj = 0. By using condition ([7]) in 
equation set ([5]), we obtain 



x^{t) =F{x,{t),m,) + aY,AtjH{xj{t),pjy 

j 

crSa,H{xi{t),p,), 



(10) 



i = l,2,...,iV, where we have introduced the matrix A' 
defined by 



A ■ ■ 

Au - Sai 



if j 7^ h 
if j = i. 



(11) 



By construction, the matrix A' = {A'^j} is such that the 
sums of its rows are constant and equal to Sa. We note 
that by setting to zero all the mismatches Smi, Spi, and 
Sai in (jlOp . a synchronized solution exists for the set of 
equations in pi7| of the form 



Xs = F^Xg, to) + crSaH{xs,p). 



(12) 
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Extended master stability function. We introduce 
the average trajectory x{t) — N^^ J^k ^k{t) that satisfies 
the following average dynamics 

k k,j 

a'^5akH{xk{t),pk) ■ 
k 

(13) 

Since the quantities 5ai, Spt, and dnii are small, we expect 
the individual trajectories Xi{t) to be close to the average 
trajectory x{t), that is, \\xi{t) — x{t)\\ < k* for all times 
and some small k* > 0. We define the variation with 
respect to the average trajectories as Sxi{t) = {xi{t) — 
x{t)). By expanding both (jB]) and to first order about 
{x{t),'m,p), we obtain 

Sxi{t) =DFx{x{t),fh)5xi{t) + DF,n{x{t),fn)Smi+ 

aDH,{x{t),p) - bj)5x,{t) + 

' (14) 
aDHp{x{t),p)Y,iK,~bj)6p,+ 

3 

aH{x{t),p)Sai, 

i = 1,...N, where bj = AJ^^-, that is, bj repre- 

sents the sum of the entries over column j of the matrix 
A' divided by N, for j — 1, ...,N. We have indicated with 
DFx and DHx the partial derivative of the functions F 
and H with respect to x, with DFm the partial deriva- 
tive of the function F with respect to to, and with DHp 
the partial derivative of the function H with respect to p. 
To obtain ([13]), we have used the properties J^j ^^j = 0: 

J2j = 0: = 0> ^Pj = 0' b] = and 

we have discarded second order terms in all the variations. 

We define £[,£'2, ^'n as the eigenvalues of the matrix 
A' . We note that since the row-sums of the matrix A' are 
equal to Sa, the matrix A' has one eigenvalue £[ — 5a, with 
associated right eigenvector Ijv- Now, we consider the 
matrix A = {Aij}, where Aij = {A'^^ — bj), and we look 

for the solutions of the eigenvalue equation Avi — XiVi. 
We observe that the matrix A has the property that both 
the sums of its rows and its columns are equal zero. Thus 
■Si = 1 AT is still a right eigenvector for the matrix A, with 
associated eigenvalue Ai = 0. Moreover, = Ijv is also 
the left eigenvector of the matrix A, associated with the 
eigenvalue 0. The remaining eigenvalues of the matrix 
A are = £'^ for i = 2, ...,N [JT]. In other words, the 
matrices A' and A have the same spectrum except for the 
eigenvalue associated with the right eigenvector vi — Ijv- 
As discussed in what follows, the eigenvalues A2,...,AAr 
control the stability of the nearly-synchronous solution. 



Equations (ITil) can be rewritten as 

dX{t) :^[In ® DF^{x{t),m) + aA® DH^{x{t),p)]5X{t) + 
[In ® DF^{x{t),ffi)]SM+ 
a[A® DHp{x{t),p)]5P+ 
a[lN®H{x{t),p)]6A, 

(15) 

where 5X{t) = [5xi{tY ,5x2{tY , ...,5xN{tYV , 6M = 
[5mi,(5TO2, ...,5mAr]'^, 6P = [Spi, Sp2, Spn]'^ , SA = 
[Sai,Sa2, ■■■,SaN]'^ , and the symbol (8> indicates direct 
product or Kronecker product. 

Following and assuming that the matrix A is 

diagonalizable, we write, A — VAW, where A — 
diag(Ai, A2, Ajv), V is a matrix whose columns are the 
right eigenvectors of the matrix A, and W = V~^. Pre- 
multiplying (IT5|) by ® /„, we obtain 

Q{t) =[/jv ® DF,{x{t),m) + (tA «) DH,{x{t),p)]Q{t)+ 
[W (g) DFmix{t),fh)]SM+ 
a[X^W ® DHp{x{t),p)]5P+ 
a[W (E) H{x{t),p)]dA, 

(16) 

where Q{t) = (W ^ In)SX{t). We note that both matrices 
In and A in the homogeneous part of Eq. are diagonal 
matrices. Thus equation ([TB|) can be decomposed into N 
blocks of the form 

q,{t) ^DF,{x{t),m) + aX,DH,{x{t),p)]q,(t) + 

Wij6mjDF„i{x{t),m) + 

akJ2w,j5pjDHp{x{t),p)+ (17) 

3 

VFy(5ajiJ(S(t)), 

3 

i = 1,...,N. We comment that the homogeneous part of 
each block in (IT71) is independent of the other blocks. For 
i = I, the variational equation (|17l) yields qi{t) — since 
EiLi ^Xi{t) — and \n is a left eigenvector. Thus we note 
that the component of the evolution along the direction 
xi = a;2 = ... = a;jv is not affected by the mismatches Snii, 
dpi, and dui. Stability of the nearly-synchronized solution 
is controlled by perturbations in the remaining directions, 
q2, ...,qN- Following [20[[3T] . it is possible to associate the 
following parametric equation to Eq. (|17p 

z{t) ^[DF,{x{t),ffi) + ujDH,{x{t),p)]z(t)+ 

eDF,n{x{t),^)+(:DHp{x{t),p)+r^H{x{t)), 

which corresponds to equation set ([TT]) upon setting z = 
qi, uj = aXi, e = J^j^ij^'^j^ C = crXiJ^jWijdpj, and 
V = <^J2j WijSttj, for i = 2, N. 

In order to assess the linear stability of the nearly- 
synchronous solution, Eq. (IT51) needs to be tested for 
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the set of eigenvalues A2, Aat. If the Lyapunov expo- 
nents associated with the homogeneous part of Eq. ([T^ , 
i = 2, ...,N are negative, the nearly-synchronous solution 
is stable. In this case, the forcing terms on the right hand 
side of Eq. p^ . i ~ 2, ...,N, can be considered as inputs 
to a stable system. It is then possible to associate an ex- 
tended master stability function e, C, 77), defined as 

lim^^oo y^r-i/J" \\z{t)\\'^dt to Eq. (UHl), which yields the 
asymptotic norm of the time average of z as a function 
of the tuple {u!,e,C,v)- However, stability of the nearly- 
synchronous solution depends on the homogeneous part 
of (|18p . that is, it depends on w, while it is independent 
of e, C, and 77. We note that for Sai — 0, 6mi = 0, and 
6pi = with i = l,...,iV, the parametric equation 
reduces to (O, which corresponds to the ideal case where 
all the parameters are equal to their nominal values. 

Moreover, following [3T], in the case that the master 
stability function is asymptotically bounded and uj is fixed, 
we have that Ai{uj,eXiV) scales linearly with respect to 
e, and rj, that is, 

Miu},eX,v) - c,iuj)\e\ + c^iuj)\C\ + c,j{uj)\r]\, (19) 

where the coefficients Ce-,Cri, and are functions of uj. 

We comment that the extended master stability func- 
tion depends on the eigenvalues of the perturbed matrix 
A' and not on those of the nominal matrix yl^OM^ rpj^g 
matrix A' can be considered a perturbed version of the 
nominal matrix A' = + A, where the per- 

turbation matrix A = {A^j} — {SAij — S'^^{J2j SAij — a)} 
and (5*^ indicates the Kronecker delta, i,j — 1, ...,N. Note 
the sums of the rows of A are equal to 5a. The eigenval- 
ues of the perturbed matrix A' can be computed from the 
spectral properties of A^^^^ . By using classical pertur- 
bation theory [33) . and assuming that the eigenvalues of 
the matrix A'^'^^^ are all distinct, we find 

A,^£fOM + ^^, ^^2,...,N, (20) 

Vi 

where Wi and Vi are the left and right eigenvectors asso- 
ciated with the eigenvalues matrix A^'~'^^ , 
respectively. Equation ([^0]) shows that the deviations of 
the relevant eigenvalues from their nominal values are 
on the same order of the perturbations A^ on the cou- 
plings. We also comment that Eq. (PH)) predicts that 
£[ ~ (i&f Aljv)/(«)f Iw) = a, since £[ = a by construc- 
tion. Similar arguments can be used to estimate the left 
eigenvectors of A from the spectral properties of A^^^^ . 

The main result of our analysis is that stability of the 
nearly-synchronous evolution for the system ([BJ can be as- 
sessed by using a master stability function, which depends 
on the eigenvalues of an appropriately modified coupling 
matrix A' . Though in a practical situation it is not fea- 
sible to exactly calculate these eigenvalues, for small de- 
viations of the couplings from their nominal values they 
differ from their nominal values ^^om small quantity 



of the same order of the A. Moreover, the mismatches in 
the individual functions F and H, along with the devia- 
tions in the row-sums of the coupling matrix A, introduce 
forcing terms in Eq. ^TE\\ through the coefficients e, C, and 
rj. Such forcing terms maintain the network in a state of 
approximate synchronization. 

Following |31| , in case the matrix A has an orthonormal 
basis of eigenvectors, that is, it is symmetric, we can write 

r N N 

E= hm T"i / V||fa,(i)fdt = Vx'(w^,e„C^,^7^). 

. 

We note that iJ is a quantity of physical interest, as it rep- 
resents the time average sum, over all the coupled systems 
of the distances \\Sxi{t)\\ from the average trajectory x{t). 
One of the advantages of this approach is that, by comput- 
ing the master stability function once, E can be estimated 
for any network topology that approximately satisfies the 
constant-row-sum condition. 

As pointed out in |31) . a complication with this ap- 
proach is that Eq. depends on x{t), which is an av- 
eraged trajectory over all the systems in the network. In 
a large network, calculating x{t) may be computationally 
expensive, as it requires full integration of N individual 
systems, see Eq. ([T^ . However, for practical purposes, 
x{t) in ([T^ can be replaced by the individual dynam- 
ics Xs(t) in (jl2l) . which depends explicitly on ffi,p, and 
6a. We comment that, unless precise knowledge of the 
characteristics of all the individual units and of their cou- 
plings is available, it is difficult to exactly compute ?7i,p, 
and Sa. Nevertheless, a priori knowledge on the statistical 
properties of the coupled systems can be used to infer the 
average parameters. For example, if mi, pi and Oi, with 
i = 1, N are taken as independent and identically dis- 
tributed random variables, drawn from distributions hav- 
ing mean corresponding to their nominal values, and finite 
variance, the central limit theorem states that rfijp, and 
da approach their nominal values as the number of nodes 
increases. 

Numerical simulation. We use the algorithm in [34] 
to generate a scale-free network of = 100 nodes with 
average degree equal to 30 and exponent of the power- 
law degree distribution equal to 3. For each pair of nodes 
i,j = 1, N, j ^ i, A^o*^ = Aj^o*^ = 1 if nodes i and 
7 are connected; otherwise, A^°*^ = = 0. We set 

Af(OM ^ which guarantees that the row- 

sums of the matrix A'^^^ are equal to zero. Moreover, 
as the matrix ^^om gymmetric, it is diagonalizable, 
its eigenvalues are real, and the eigenvectors can be taken 
to be orthonormal. We find that £2 = 12.7018 and In = 
86.0531, where we set £1 < £2, < £n- 

We consider that Aij = Afj'^^ {1 + <;aPij), for i,j = 
1, N, where pij = pji is a random number drawn from a 
standard normal distribution and <ia is a scalar. Upon this 
selection, the matrix A' in (jll[) is symmetric; the matrix 
A is also symmetric, since bj — Sa/N for j = 1, A^. For 
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= 10-4, we obtain A2 = 12.7007 and Aat = 86.0529. 
This is in agreement with Eq. ([^0]). as we find that jA^ — 
ii\ is on average on the same order of magnitude of the 
deviations on the couphngt0. 

We perform a numerical experiment for a set of nomi- 
nally identical Rossler oscillators that are affected by mis- 
matches in both their dynamics and output functions and 
are coupled by the scale free network, described by the 
matrix A. In this case, the equations of motion are 



iiiii) = - Xi2{t) - Xi3{t) + a^Aij{xji{t) +Pj), 

j 

Xt2{t) =Xii{t) +miXi2{t), 

iait) =0.2+(a;,i(t)-7)x,3(t). 



(22) 



i = 1, N, where the state vector of oscillator i is Xi — 
[xii,Xi2,Xi3]'^ . The parameters pj are random numbers 
drawn from a Gaussian distribution with mean zero and 
standard deviation c^p, and the parameters rrii are random 
numbers drawn from a Gaussian distribution with mean 
value equal to 0.2 and standard deviation <^„i. 

In Fig. 1, we plot the error measure E, defined in (HI]), 
versus the coupling strength a. Simulations are run for 
a total time duration T = 3000, which is considerably 
larger than the typical time scale of an oscillation for an 
uncoupled Rossler oscillators, that is 27r; time averages are 
taken over the time interval [2700,3000]. 

From the direct numerical integration of Eqs. and 
(fTS]) with TO = 0.2, p = 0, (5a = 0, and e = C = ?7 = 0, 
we find that the master stability function converges to 
zero in the range 0.143 ^ a; ^ 4.40, which for our choice 
of the matrix A, corresponds to stability in the range 
0.0113 <(T < 0.0511. This range is delimited by the ver- 
tical dashed lines in Fig. 1, which shows good agreement 
with our computations of the full nonlinear system ((22|) . 
Figure 1 illustrates that the range of stability is affected 
neither by the presence of small deviations from the nomi- 
nal couplings nor from small mismatches in the individual 
oscillators' parameters. This is because the eigenvalues Xi 
are indistinguishable from the eigenvalues £i, for i = 2 or 
TV to the degree of accuracy of the simulation shown in 
the figure. However, for a inside the range of stability, the 
value attained by E depends on the values of Sui, 5mi, 
and 5pi. Figure 2 shows c^, cq, and versus w. With this 
information, Eq. (|19l) provides an estimate of the mas- 
ter stability function for any tuple (w, e, C, ??)■ We use Eq. 
(fT9|) along with the data plotted in Fig. 2 to calculate 
the master stability function M.. This is shown for com- 
parison in Fig. 1, where the symbols x {+) are used to 

plot Ya=.2 J^i.'^ii ^iX'i, "Hif for = 10"'' and <r„i = <rp = 
{<,a = 10-4 and — — ^ X 10^4) _ Poorer agreement is 
observed for values of a slightly above the lower threshold 
for stability of 0.0113 (not shown), which corresponds to 



■^We have also performed numerical experiments for pij ^ pji and 
we have found that, for small values of fa, the eigenvalues A^'s are 
still real. 
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Fig. 1: Triangles, diamonds, and squares represent the error 
measure E versus the coupling strength a. Triangles are used 
for the case in which <;a = ?m = ?p = 0. Diamonds are used for 
the case in which <;„ — and <;m = ?p = 0. Squares are used 
for the case in which ~ 10-^ and w = ?p = 5 x IQ-*. The 
vertical dashed lines delimit the range of stability predicted 
by the master stability function. The symbols x (-f ) refer to 

I]^2-'^('^"<^i'Ci.'70^ for ?a = 10~* and w = ?p = (?a = 
IQ-'' and w = ?p = 5 x IQ-*), computed using Eq. (|19p . 



a so-called bubbling region, as further discusses belowo. 

Conclusions. The master stability function analysis 
[19l[20] provides a necessary and sufficient condition for 
linear stability of the synchronous solution for an arbitrary 
network of coupled identical systems. An extension of this 
approach for networks of groups, where the dynamics of 
nodes within a group are the same but are different for 
nodes in distinct groups, is proposed in [35]. In addition, 
a master stability function for networks in which each unit 
independently implements an adaptive strategy to main- 
tain synchronization is presented in [36 . The analysis of 
nearly identical coupled dynamical systems is considered 
in [5U1I5T| . For this case, which is of practical relevance in 
experimental instances of synchronization and in biolog- 
ical systems, it is shown that a master stability function 
approach is applicable [3Tj. 

In this letter, we have proposed a sensitivity analysis to 
address synchronization in the presence of a broad range 
of deviations from nominal conditions. In particular, we 
have taken into consideration simultaneous small devia- 
tions in the dynamics of individual units, the output func- 
tions of the individual units, and the coupling among the 
systems. We have shown that the master stability function 
formalism can be extended to this general scenario and 
that stability of the nearly-synchronous evolution depends 
on the eigenvalues of an appropriately modified coupling 



^ Numerical experiments performed by replacing Snii, Spj, 
and SAij's with random numbers from the same distributions and 
Ai, Ajv and W with the eigenvalues and eigenvectors of the orig- 
inal matrix A'^'-'^^ show good agreement with the results in Fig. 
1. 



p-5 



F. Sorrentino and M. Porfiri 



400 r 
200 - 




0.5 1 1.5 2 2.5 3 3.5 4 




0.5 1 1.5 2 2.5 3 3.5 4 



400 n 
200 - 
0- 




0.5 1 1.5 2 2.5 3 3.5 4 



Fig. 2: Ce, cq, and c,, versus ui. 

matrix. Our analysis is motivated by inherent practical 
challenges in implementing ideal conditions in experimen- 
tal analysis of synchronization. For example, our approach 
can be directly applied to synchronization of nearly iden- 
tical units whose interconnections yield to approximately 
zero-row-sum coupling matrix. In this case, the proposed 
master stability function can be used to estimate the con- 
ditions under which the nearly-synchronous evolution is 
stable and in case of stability, the approach can be used 
to quantify the overall synchronization error. 

Noise or small mismatches in the parameters of the in- 
dividual systems can be responsible for the onset of hub- 
hling [5ni[551[57] , that is, rare intermittent large deviations 
from synchronization. We expect bubbling also to arise 
in the case of approximate satisfaction of the zero-row- 
sum condition; in this case, the master stability function, 
introduced in this letter, can be used to identify stable, 
unstable, and bubbling regions in the relevant parameter 
space, see for example |36) . 

* * * 
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